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ABSTRACT 



This paper introduces methods tailored especially for problems 
whose solution behaves like e'“, where A is complex. The shallow 
water equations with topography admit such solution. 

This paper complements the results of Pratt and others on 
exponential-fitted methods and those of Gautschi , Neta, van der 
Houwen and others on trigonometrically-fitted methods. 

1 . Introduction 

In this paper we consider linear multistep methods 
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3 - y 1 - 

A,' n+l- V, 



= h 



k 

i = 0 






1 , n 



k-1 



( 1 ) 



for integrating the initial value problem 



y'(x) = f(x,y(x)), v(Xq) = y^ . 



( 2 ) 



This linear .multistep method is characterized by the polynomials 



P ( ■; ) = ) a , c 
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The main assumption of this paper is that it is a priori known 
that the solution is appro.x ima te ly of the forri 
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where ” '^j j ' frequencies w^ are in a given interval 

(w^,WuJ . 

The special case where = jw^ with w^ given was considered 
first by Gautschi [8]. His approach was the following. Let: 



d(z) = p(e^) - zo(e^) 



( 5 : 



then the local truncation error of (1) is given by Lambert [11] 
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Insertina (4) in (6) vields 
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The coefficients b^ are chosen in such a way that 

t(ihjwQ) =0,j=0,l,...,q, (8) 

for the largest value of q possible. q is then called the trigo- 
nometric order of the method. Gautschi has chosen a. such that the 

1 

methods are of Adams and Stormer type. However, these methods are 
sensitive, to changes in the frequency w’q . Neta and Ford [13] 
developed Nystrom and generalized Milne-Simpson type methods. These 
methods showed less sensitivity to perturbation in w^ but require 
the eigenvalues of the Jacobian to be purely imaainary. Neta [14] 
has developed families of backward differentiation methods that 
overcome the above-mentioned restriction. Salzer [17] has developed 
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predictor-corrector methods based on trigonometric polynomials. 
See also Steifel and Bettis [18] and Bettis [3] . Van der Houwen 
and Sommeijer [10] have developed an alternative approach. The 
conditions (8) were replaced by 



( 0 ) = 0 , 



(9) 



d (ihX = 0 , j = 1,2, . . . ,q , 



where the are appropriately chosen points in the interval 

[w^,Wu] . 

An advantage of this so-called m.inimax approach over the fitting 
approach is the increased accuracy in cases where no accurate esti- 
mate of Wq is available or when the frequency is varying in time. 

The other special case considered in the literature is where 
Aj = i’- j • Probably the first article on the subject is due to Brock 
and .Murray [5]. They discuss the use of exponential sums in the 
integration of a system of first order ordinary differential equatio.ns 
Dennis [7] also suggested special m.ethods for problems whose solution 
is exponential. He suggested a transformation of variables. .More 
recently, Carroll [6] has developed exponentially fitted one-step 
methods for the scalar Riccati equation. For the general first 
order system of equations, Pratt [16] suggests methods based on the 
three parameter exponential function 



I (x) = A + Be 



zx 



(lo: 



H 



The parameters A, B are given in terms of values of y and f. 

Several possibilities for z are given based on results of Brandon 
[2] and Babcock et al . [1] . 

Lyche [12] analyzes multistep methods which exactly integrate 
m '^n^ 

the set {x e }, where w^ is real or imaginary. 

In this article we developed various methods fitting exponentials 
and methods obtained via the minimax approach. 

2 . Construction of Methods 
2 . 1 Fitting Methods 

In this subsection we discuss various fitting methods. To this 
end, we separate i(ihjX) = 0, j = l,2,..,,q, into real and imaginary 
parts. This yields the following equations relating the coefficients 
a, ,b, , 



''^cos jv(k-i) 

v=0 ■' 



I a^e^"^^" ^'^sin jv(k-i) 
£ = 0 



J b„e^’’^^^ ^ ^ [ ju cos jv (k-il) 

:; = 0 "" 

- jv sinjv(k-£)] = 0 , 

k . 

I b„e^^ [ jy sinjv(k-i) 

x = 0 

+ jv cosjv(k-f)] = 0 , 



where >. = w -h iip , v = ~hC', v = hw , j = l,2,...,q. 

For explicit methods, b^ = 0. For Adams type methods a^ = 1, a^^ = -1 
a^ = 0 for i = 2,...,k. For Nystrom or generalized Milne-Simpson 
methods a^ = 1, a 2 = -1 and other a^ = 0. 
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k = 1 Implicit 
Adams 



^0 = 



ve ^ + 



sinv - 



cos V 



2 



sin V 



b 



1 



ve^ 



- U sinv - 
sin 



V cos V 



V 



l — C'C) c: \» 

For i' = 0, the coefficients become = : = b 

with Gautschi [8] if the coefficients are expanded in 
with respect to v. 



k = 2 Explicit 
Adams 



^1 = 



(u sin2v - V cos2v)e“ + v cosv - i sinv 

2 2 

(u +■•- ) sin V 



^2 = 



(v cosv - v sinv)e “ - ve‘ 
2 

(i +v‘") sin V 



Nv Strom 



b, = 



(y sin2v - v cos2v)e"' ve 



2 2 

(u +v ) sin V 



^2 = 



(v cosv - y sinv)e^"' + ve 
2 2 

(u +v ) Sin V 



„ , , .. -.sinv. 

For 1^=0, the coefficients become b^^ = 2 — ~ — , = 

agree with Neta [14] . 



( 12 ) 

, which agree 
Taylor series 



( 13 ) 



( 14 ) 
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k = 2 Implicit 



In this case, one obtains a one-parameter family of (Adams, 
generalized Milne-Simpson ) methods of trigonometric order 1. The 
free parameter can be used to increase the algebraic order of the 
method as in [13] . 

Backward Differentiation 

e^^ cos 2v + a^e^ cos v + a^ ~ bQe^^(u cos2v - v sin2v) = 0 , 

e^"^ sin 2v + ^ ~ sin2v + v cos2v) = 0 , (15) 

1 + a^ + ^2 = 0 . 

This system can be solved by MACSYMA (Project MAC ' s SYmbolic MAnipu- 
lation system written in LISP and used for performing symbolic as 
well as numerical m.athematical manipulation [4]) or by REDUCE [9]. 
The solution is 

2u . _ _ 

ve - u sin2v -v cos2v 

^1 ~ ” ' 

-e" (v cosv + u sinv) + (j sin2v + v cos 2v 

(16) 

a2 = -1 - 3i , 

2v . u • T 

-e sin V + e sin2v - sin v 

e^^^(v cosv + sinv) + (^j sin2v + v cos2v) 

For = 0 , the coefficients agree with those given in [14] . 
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k = 3 Explicit 



Again here, one obtains a one-parameter family of methods of 
trigonometric order 1. In order to get methods of trigonometric 
order 2, one has to construct a 3 step implicit method of Adams or 
generalized Milne-Simpson type. In order to increase the trigo- 
nometric order without going to a higher step number, one can con- 
struct linear multistep methods for which the coefficients a^ are 
also functions of , w. Some examples are given in the next 
subsection . 

2 . 2 Generalized Fitting Methods 

In this section, we construct some linear mulstistep methods 
of the form 

k k 

V 

Since a, are functions of A one has more free parameters for his 
disposal which can be used to obtain higher trigonometric order 
methods with relatively lower step number. 

k = 2 Implicit 

In this case, one has to solve the following linear system of 
five equations for the parameters a^ , a 2 , bQ, b^ , b 2 to obtain a 
method of trigonometric order 2. 



a. (,\)y = h [ b, (A) 

0 '■ "■I--- i=0 - 



n^l-£ 



k > 1, n > k-1 



(1 
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1 + ^ ^2 ~ ' 

cos2v cos v "*'^2 cos2v -v sin2v) 

-bj^e^(vi cosv -V sinv) - = 0 , 

sin2v +3^0^^ sinv -bQe^'^(y sin2v +v cos2v) 

-bj^e'^(u sinv +v cosv) -vb 2 = 0 , (18) 



4u , ^ 2u 

e cos 4 V + a j^e 


cos 2 v+a 2 ~ b^e^^ 


(2y 


cos 4v 


-2v 


sin 4 V ) 




-b^e 


(2y 


cos 2 V 


-2v 


sin 2 V ) -2yb 


4u . . 2y 

e sin4v +a^^e 


sin2v -b^e'^"' 


(2u 


sin 4v 


+ 2v 


cos 4 V ) 






(2u 


sin 2v 


+ 2v 


cos 2 V ) -2vb 



The system was solved by REDUCE [9]. The expressions for the 
coefficients are complicated but REDUCE produces an output in the 
form of Fortran statements that can be incorporated in a computer 
program, for numerical experiments with such a method. 

2 . 3 Minimax Methods 

In this section we discuss minimax methods, i.e., methods 
obtained by satisfying conditions (9). These conditions can be 
written in terms of a„, b„ as follows: 



I a.e^’"' ^cos (k-£) V ^ 



k 

I 

J?. = 0 



= " b.e^^-')’^^'^-v(^)sin(k-;:)v(^^u(^)cos(k-Uv(^^}, 

i_ = 0 " 



( 19 ) 



j ~ , 



9 



ik-l) u 



( 3 ) 



I a„e''' sin(k-£)v 



(j) 



£ = 0 



^ / 1 
= I b.e 
£=0 



-M- (3) 



(j) 



(j) 



{v-^cos(k->. )v-^+U-'sin(k-j^) 



(3) . 



( 20 ; 



j = 1,2, 



where r'^* = h£ '3> , v'i> = hw'^) . 

are the zeros of the function c(ih"') such that it has 
a small maximum norm in the rectangle ^ £ ■* ^ -u’ 

To obtain the best approximation in this case is certainly 
not easy; but we will assume that one can write ^(ih'^ ^^^) = 
as a product of 2 one variable functions. Thus 
and w^^^ can be taken as Chebyshev ' s points on the corres- 
ponding interval, i.e. 



(j) 




- Y cos 



2 j-1 

2q 



w 



(j) 



1 , 

2 (w^, 



+ w^) 



■^(w^ -w^)cos 



2n-l 

2q 



j ~ 1^2, ...^q 



( 21 , 



For this choice of points, one can evaluate the coefficients 
a „ , b^. by solving 

X. Aw 



: (0) = 0 , 

; ( y ^ ^ h w ^ ^ M = 0 , j = 1 , 2 , . . . , q . 

. • 2 

We call such methods product minimax (pM ) . 



( 22 ) 
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The number of free parameters for implicit methods 2k+l and 
the number of equations is 2q+l, thus, the trigonometric order 
q is equal to the step number k. 

k = 1 Implicit 

In this case 






and the system of equations can be solved for b^, b^^ . This 
yields the coefficients given by (12) where 

M = , V = hw^^^ . (24) 

Thus, the product minimax method would suggest using the center 

of the rectangle I , '|^^] x [w^ , w^] as Aq. To obtain a product 

minimax method of trigonometric order 2, one has to solve a 

system of 5 equations similar to equation (18) with the unknowns 

b^, bj^, b 2 , , a 2 • The difference is that in the last 2 equa- 

( 2 ) ( 2 ) 

tions one should replace 2y by y and 2v bv v . In the 

second and third equations of (18) , the y, v, should be replaced 
by y^^^, respectively. The resulting system can be solved 

by MACSYMA [4] or REDUCE [9]. 

In the next section we implement two methods of trigonometric 
order 1 and 2 and see how the product minimax methods compare 
with fitting methods. 
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3. 



Numerical Example 



Both systems (18) and (19) -(21) were solved by REDUCE which 
produced a FORTRAN subroutine for the evaluation of the coeffi- 
cients. This subroutine is called only once during the integration. 

The methods were compared for the solution of the initial 
value problem 



z - j i(l +i)z = 0 , ® 1 

z(0) = 1 



(25) 



whose exact solution 



iy(l+i) t 
z(t) = e ^ 



( 26 ) 



thus 



A = ^ ~ 2 ' " 2 ’ 

In order to avoid complex arithmetic, we rewrite the differential 
equation as a system of equations for the real and imaginary 
part of z = u + iv. 

u + ^(u+v) =0 , 

0 ^ t £ 4 

v--^(u-v)=0, (28) 

u (0) = 1 , 

v(0) = 0 . 
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The system is solved by fitting methods of trigonometric order 
1 and 2 with h = .01 and various values of \p and w. In Table 
1 we list the Euclidean norm of the error at t = 4 . It is 
clear that the method is not sensitive to perturbations in the 
values of ip and w. 





Aw 


error 


first order 


second order 


0 


0 


. 3678 (-12) 


.4433 (-12) 


0 


.1 


.4482 (-7) 


.3508 (-11) 


.1 


0 


.4346 (-7) 


.2544 (-11) 


.1 


.1 


. 6342 (-7) 


.4421 (-11) 


0 


. 2 


. 9241 (-7) 


.8126 (-11) 


.2 


0 


. 8706 (-7) 


.5095 (-11) 



Table 1 

Using the product minimax methods of trigonometric order 
1 and 2 with h = .01 and various squares centered at = - tj/2, 
w = tt/2, the error is much laraer but again is insensitive 
to small perturbations in the length of the sides of the 
squares. In Table 2 we list the Euclidean norm of the error 
at t = 4 . 
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length of side 


r— r 

error 


first order 


second order 

■ 1 


.4 


. 3678 (-12) 


.1469 (-6) 


. 8 


. 3678 (-12) 


.2343 (-5) 


1.2 


. 3678 (-12) 


. 1186 (-4) 


1.6 


. 3678 (-12) 


. 3728 (-4) 


i 2 

! 


. 3678 (-12) 

1 


.9001 (-4) 

1 



Table 2 

Note that the perturbations in the product ninimax methods are 

larger than those allowed in the fitting methods. It is 

possible that the larger errors in the product minimax methods 

are due to the assum.ption that ■: can be v.’ritten as a product 

of 2 one variable functions. Also note that for the first 

2 

order method, one always gets a good result since PM always 
uses the center of the square. 
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APPENDIX 



Here we show that the shallow water equations witn 
topography have a solution of the form e^"", where A is complex. 
This system of equations consists of three equations with 



9u 



3 V 



variable 


2S U , 


V 


and (p . 


The equations 


9u 


3u 


fv 


3d; 






9y 


1 u V 

dx 


^ 0 , 


9v ^ 


? V 


fu 


3 i'*' 






dy 


+ = 
dy 


= 0 , 



(A.i; 



(A. 2 : 






(A. 3) 



where - gh is the geopotential height (h = height of free 
surface) , is the bottom topography (assumed to be independent 
or time) , u, v are the components of the wind velocity in the 
X, y direction, respectively, and f is the Coriolis parameter. 
Linearizing the equations by letting 

u=U+u', v = V + v', 0 = 

where U, V are the constant mean flow and b is independent of 

time. Assuming that U, V are related to $ via the geostrophic 
relations 



U 



f Sy ' ^ ~ I ^ 



(A. 4) 



one obtains the linear system (after dropping the primes) : 
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+ U + V - fv ^ = 0 , 

3 t dx dy dX 



(A. 5) 



+ U + V + fu + = 0 , 

3t dX dy ay 



(A. 6) 



|i + u + V H + iiHXl + iiTLl = u V 

3t 9x 3y dX dy oX oy 



dC- 



B 



(A. 7) 



where y = 4> - 4< 



B' 



If the flow is assumed to be along the topography as in 
[15] , then the right hand side of (A. 7) is zero. In such a 
case, one can write the solution in the form 



u = u e 
o 



i ( fx + ^y - ct ) 



V = V e 
o 



i(fx + rv - "■( 



(A. 8) 



i (:.x +-y - ct) 



where 



r 



ip , 



(A. 9) 



n = V - i 6 . 



In order for (A. 8) to be a solutio.n for (A.5)-(A.7), one must 
have 



A -- — o + + gV 



(A. 10) 



satisfying 
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iA [ (iA ) - in (iny + 



- 

3yJ 



f [-iA f - ir (inY + 



iX) 1 

9y^ ■' 



+ 



(iCY 



+ 



lx 

X 



(-inf 



+ CA) 



0 . 



(A. 11) 



The real and imaginary parts of 



A^ = -a + yU + vV 

A^ = -pu - ev , 



(A. 12) 



satisfy the following system of equations (after dropping 
nonlinear terms in A) 



2 2 2 
[f^ +Y(l +n )] 






3 Y > x: , 3 Y 

n 7T-7) = f (n 

oy dx 



- c it) 

dy 



A„ (< 



lx + . IX) 

3x ^ 3v' 



- > . [: 



+ Y(C^ 



1 

n ) J 



= 0 . 



(A. 13: 



In general, / is complex and, thus, the shallow water equations 
have a solution in the class of problems to be discussed here. 
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